;________________________________________________________________
; --> FOR MIXING LENGTH:
;________________________________________________________________
;
;  The reference comparison to be made with MLT is given by:
;
;          DT/T ~ <w^2>/cs^2 ~ [F_conv/(rho_strat*cs^3)]^2/3 
;
;  The relative temperature fluctuations are given by:
;
;                      DT/T = <Th'>/(<Th'>+theta_strat)
;
;  the speed of sound is given by the usual definition:
;
;                      cs=sqrt(dP/d(rho))
;
;  which for this case is given by the polytropic atmosphere
;  and so we get:
;  
;                      cs=sqrt((1+1/n)*(P_strat)/(rho_strat))
;________________________________________________________________

;________________________________________________________________
; UPDATE THE POLYTROPIC INDEX PARAMETERS !!
;________________________________________________________________

if m1 eq !Null then begin
 print, 'WARNING !!! '
 print, 'Set the parameters of the polytropic index profile used: '
 print, '(m1 , m2 , m3 , xtac , xtop , d) ... '
endif else begin
 mr=m1-(m1-m2)*0.5*(1+erf((z-xtac)/(d)))-(m2-m3)*0.5*(1+erf((z-xtop)/(d)))

;________________________________________________________________
; COMPUTE THE SPEED OF SOUND -> polytropic profile
;________________________________________________________________

 cs=sqrt((1+1/mr)*rg*temp_strat)

 print, 'INSTANTANEOUS PROFILES TAKEN AT topslice'
 print, 'current value of topslice = ', topslice
 ;________________________________________________________________
 ;Reference dimesionless mixing length quantities taken at topslice
 ;________________________________________________________________

 g1=make_array(z_num)
 g2=make_array(z_num)
 g3=make_array(z_num)
 
 for jj=0,z_num-1 do begin
  g1(jj) = the_aver(jj,topslice)/(the_aver(jj,topslice)+theta_strat(jj))
  g2(jj) = w2_aver(jj,topslice)/((cs(jj))^2.)
  g3(jj) = fconv(jj,topslice)/(rho_strat(jj)*(cs(jj))^3.)
 endfor

 ;________________________________________________________________
 ;TIME AVERAGED mixing length quantities
 ;________________________________________________________________
 
 g1_prom=make_array(z_num)
 g2_prom=make_array(z_num)
 g3_prom=make_array(z_num)
  
 for jj=0,z_num-1 do begin
  g1_prom(jj) = the_prom(jj)/(the_prom(jj)+theta_strat(jj))
  g2_prom(jj) = w2_prom(jj)/((cs(jj))^2.)
  g3_prom(jj) = fconv_prom(jj)/(rho_strat(jj)*(cs(jj))^3.)
 endfor

endelse
end
